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Abstract 



o 

^^ We consider a weakly nonlinear solution of the Cauchy problem for the regularised Boussinesq 

^> equation, which constitutes an extension of the classical d'Alembert's formula for the linear wave 

^ equation. The solution is given by a simple and explicit formula, expressed in terms of two special 

^ functions solving the initial- value problems for two Korteweg-de Vries equations. We test the formula 

^^ by considering several examples with 'exactly solvable initial conditions' and their perturbations. 

Explicit analytical solutions are compared with the results of direct numerical simulations. 

^^ Keywords: Boussinesq equation; Initial-value problem; Asymptotic multiple-scale expansion; 

'; Solitons; Averaging. 

a 

1 Introduction 

(N 

I The Boussinesq equation and its generalisations, since the original derivation in the context of fluids 

lO [H] and reappearance in connection with the famous Fermi-Pasta-Ulam problem |141l42j. have recently 

r emerged in a vast variety of problems describing nonlinear waves in solids (see, for example, \W\ [321 

[371 [361 EH I2EI US]). At the same time, considerable progress has been made in understanding the 

validity of regularised models and proofs of existence and local well-posedness of the initial-value 

PsJ problem (IVP) in the context of water waves (see [H [5l |6l [71 [31] and references therein), while the 

^~~' global well-posedness is known to be a complicated issue [23]. Some progress also has been made 

t> in the study of boundary- value problems [15j . Most relevant to our present paper are the results 

k> establishing the validity of two Korteweg-de Vries equations (KdV) as a leading order approximation 

Vh to solutions of Boussinesq- type equations, as well as some results for the higher order corrections 

^ [m dB [251 [211 [221 [381 E [Ml 131] • 

Recently, we constructed a weakly nonlinear solution of the initial- value problem for a system of 
coupled Boussinesq equations on the infinite line [27] , assuming that initial data generates sufficiently 
rapidly decaying right- and left-propagating waves, in terms of solutions of various Ostrovsky-type 
equations (see j35j for the original Ostrovsky equation). When coupling parameters are equal to 
zero, our results yield a simple formula for a weakly nonlinear solution of an IVP for the regularised 
Boussinesq equation 
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also known as the improved Boussinesq equation. The advantage of using model (II]) is that it obviates 
the short-wave instability (see [5), which is inherent in other versions of the Boussinesq equation, 
although this version is not integrable |43| [1] . 

The paper is organised as follows. In Section [2] we firstly generalise our formula obtained in 
|27| by considering the case when initial conditions are split into 0(1) and 0(e) terms, allowing one 
to construct explicit analytical solutions for a wider class of initial conditions. We use asymptotic 
multiple-scale expansions, similar to the type used in a study of oblique interaction of solitary waves 
[33], and an averaging, used for example in f3t f27], to derive asymptotically exact models which 
describe the leading order terms in the expansions. The leading order terms satisfy the IVP for 
two KdV equations [S^, integrable by the Inverse Scattering Transform (1ST) [IH (see also pi fT^)- 
We use two arbitrary functions present in the higher order terms to improve the accuracy of the 
approximate solution. The derived formula constitutes a weakly nonlinear extension of the classical 
d'Alembert's formula for the linear wave equation, and has a similar structure. To the best of our 
knowledge such a formula was not suggested in previous studies. The simple and explicit form of 
this solution allows us to construct a large class of approximate solutions corresponding to 'exactly 
solvable initial conditions' (see Section HI). The weakly nonlinear approach used in this paper can 
be applied to any form of the Boussinesq equation and coupled systems of Boussinesq equations (see 
|27| ) and, in particular, it offers an alternative to implementing the 1ST to the integrable version of 
the equation in physically relevant cases. 

In Section ^ we introduce a finite difference scheme to solve a scaled version of (II]) , in the form 

Jit ~ J XX ~ 2^J )xx + Jtixxi \^) 



where f = ef, t = y^ t, x = -^/e x. Equation (|2]) has particular solitary wave solutions of the form 

f = A sech^ (^^] , where A = 3{v'^ - 1), A = -y^^- (3) 

V ^ / V f ^ — 1 

We analyse and discuss the schemes stability and directly measure its accuracy by comparing the 
analytical solution ^ (denoted fsoi) to corresponding numerical simulations (fnum), at times within 
the derived asymptotic models validity (i.e. from i = to some point between i ^ ta = —j= and 

l^h = ^). 

In Section [4] we explicitly derive the weakly nonlinear solution for exactly solvable initial conditions 
of the IVP in the form of right-propagating and both right- and left-propagating N-soliton solutions of 
the KdV equation. Note that these initial conditions do not correspond to the exact soliton solutions 
of the Boussinesq equation. 

In Section [5] we consider particular cases of initial conditions, namely one- and two-soliton so- 
lutions of the KdV equation, and analyse the absolute error in comparison with relevant numerical 
simulations. We define the maximum error over x of the solution at t = r as 

e\= max |/n«m(5, r) - //(x,r)|, for / = 1,2, 
—L<x<L 

where we restrict x to the finite domain 2L and fi{x,t) denotes the weakly nonlinear solution up to 
and including O(e') terms. 

In Section|6]we consider perturbations to exactly solvable initial conditions of the IVP, in particular 
for the case of right-propagating N-soliton solutions, and again analyse the error for a specific example. 
We finish with concluding remarks in Section 7. 



2 Weakly nonlinear solution 

We consider the following IVP for a single regularised Boussinesq equation on the infinite line 

1 
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f\t=o = F^{x) + eF\x) + 0{e^), ft\t=o = V\x) + eV\x) + 0[^), (4) 

(or any other asymptotically equivalent form of this equation) for the case when the initial conditions 
generate sufficiently rapidly decaying right- and left-propagating waves. Therefore to leading order 
the initial {t = 0(1)) evolution of the Cauchy data is described by the classical d'Alembert's solution 

h{x.t) = f^{x-t) + f^{x + t) for f^{x±t) = UF\x±t)± j"" V\x)dxY 

In general /q are some step like functions but in what follows we shall restrict our considerations to the 
case when these functions are sufficiently rapidly decaying at infinity (in particular, J_ V{x)dx = 0). 
To describe the subsequent {t = 0(e~^)) evolution of the Cauchy data we introduce the slow time 
variable T = et and seek the following weakly nonlinear solution in the form of asymptotic multiple 
scale expansions 

f = r{tT) + f+iv,T) + ef\tv,T) + Oie^), (5) 

where S, = x — t and rj = x + t. Substituting ([5]) into the Boussinesq equation Q we find at leading 
order the equation is satisfied, whilst at 0(e) we obtain 

-4/1^ = (2/^+r/-+/-^^)5 + (-2/+ + /+/+ + /+^), + 2/-/+ + /+/-^ + r/+. (6) 



The subsequent derivation can be performed either by integrating (rol) and requiring that /^ is non- 
secular (see, for example, [17]) or using averaging arguments (see, for example, O EZ])- In what 
follows we use the latter of the two. Namely, we average ([6]) with respect to the fast time variable t 
at constant ^ or rj; the left hand side of ([g]) at constant S, is averaged as follows 

hm - r fldt =\mi^ r^ f^drj = lim ^ [/ll'lf " = 0, 

where we assume /^ and its derivatives remain bounded (required to have a nonsecular expansion 
([5)), and similarly we get zero when averaging the same term at constant rj. Averaging entirely over 
m) at constant ^ and assuming /^ and their derivatives remain bounded and are sufficiently rapidly 
decaying at infinity for any fixed T (consistent with relevant numerical experiments), we obtain 



Averaging entirely over ([6|) at constant rj we derive a similar equation for the function /+ . Integrating 
each of these equations with respect to their respective characteristic variables, and taking into account 
the behaviour of / at infinity, yields the following two KdV equations 



r,/ J^ ~^ 9-'€C€ — U, Jrp r)J Jri r^JirnV 



/t + orf[ + 7t/£.£ = 0, /^ - -/+/+ - -/+ = 0. (7) 
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The higher order correction /^ is then obtained by substituting ([7]) into ([6]) to yield 



(8) 



Finally substituting the weakly nonlinear solution ([5| into the initial conditions of the IVP Q, we 
derive at leading order initial conditions with respect to T, for the leading order terms / 



/ \t= 



fi = \{F\x±t)± p\\x)dx^ , 



whilst at 0(e) we obtain, within the accuracy of the problem formulation, the following d'Alembert's- 
like formulae for the functions (j){^,T) and ip{r],T): 
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where 
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+ F\x), 



t=o 



fr+fr+l (/+/f - rfv + /?€ / ^^^^ - fvv J /"^^ 



V\x). (10) 



t=o 



The dependence of the functions (p and ^ on T is inherited from their dependence on the leading order 
functions /~ and /"*", however this may be neglected, at least for sufficiently small values of time. To 
construct a more accurate solution, valid for greater values of time, one needs to know higher-order 
terms in the problem formulation Q. 

3 Numerical scheme 

We next discuss the numerical scheme used to solve the Boussinesq equation Q. We implement 
a finite difference scheme derived in 0O| for a regularised Boussinesq equation or equivalently the 
scheme used in [27] for the reduction g = 6 = 'y = and c = a = f3 = 1. 

We let X € [—L, L], for finite L, and discretise the (x, t) domain into a grid with spacings Ax = h 
and At = k. The solution f{x, i) of ([2| is approximated by the solution /(z/i, j'k) (for i = 0, 1, ..., N 
and j = 0, 1, ...) of the difference scheme, denoted fij. 

Substituting central difference approximations into equation ([2| we derive the following difference 
scheme 

— fi-i,j+i + {2 + h )fi,j+i — fi+i,j+i = (k — 2)[/j_ij — 2/jj + fi+ij] + 2h fij 

+fi^lj-l - (2 + h^)fi,j-l + f^+l,J-l■ (11) 

We choose the boundary conditions such that they are sufficiently far away from the propagating 
waves. Periodic boundary conditions were also considered but, for the added complexity and compu- 
tational effort involved in solving the difference scheme, the difference in the error of the solution was 
negligible. Thus we set 



fo,j = fN,j = 0, Vj. 



Initial conditions for simulations of the scheme are chosen to coincide with initial conditions in the 
IVP Q . However for the purpose of analysis of error in this section we choose the initial conditions 
in the form of the particular Boussinesq soliton solution ([3]) and thus 

/j,o = A sech^ ( A ) ' fiA = ^ sech^ f — - — j , Wi, 

where A = 3(w^ — 1), A = 2v{v'^ — 1)^5 and we choose v ~ y^f + 1 to ensure the amplitude is 
0(e) (required to make comparisons with the weakly nonlinear solution in the next sections). The 



nine point implicit difference scheme ( 11 ), with tri-diagonal matrices of constant coefficients, is solved 
using a Thomas Algorithm (e.g., [2]). 

To examine the schemes stability we firstly linearise by setting fij = fo+fij, where /o is a constant 
such that /o > fij Vi, j. Then using a Von-Neumann stability analysis we substitute fij = G-'e'^*^ 
(where i^ = —1) into the linearised version of (11) and derive the condition 



„ 2/c2(l + /o)sin2f 

G^ - 2^G + 1 = where ^ = 1 - ^ ^ 



/i2 + 4sin2f 
For stability we require \G\ < IVO and arbitrary k, h, which is true provided |^| < 1 and thus implies 



//i2 + 4 

Hence the roots of the quadratic in G have modulus one and the linearised form of the difference scheme 
( |ll| is stable provided k < Kc- (In practice, we used a stricter condition k < \kc, to accommodate 
for the effects of nonlinearity) . It can be shown that the principal truncation error of the linearised 
scheme is 0(/i^k^ + /i^k^). 



For a given h the restriction (12) can be used to determine valid choices of k to ensure stability. 
Table 1 displays the maximum absolute error across x for different choices of n for fixed e (and 
hence fixed v) for two different space step sizes /i, and at two different times ta and th (within the 
asymptotic region of validity of the KdV models ([T])). The optimal discretisations h = k = 0.1 are 
chosen for simulations where the maximum absolute error within the time interval considered ranges 
from O(10~^) to O(10~^). This magnitude of error of the numerical solution is deemed suitable from 
previous numerical studies on Boussinesq- type equations (for example, see [20], [M!, [E], [9], [TH]). 

Max error at t = 32 Max error at t = 300 

K h=0.1 h=0.01 h=0.1 h=0.01 
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4.4693 X 10-^ 


4.5089 X 10-^ 


3.3 X 10-3 


3.3 X 10-=^ 
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8.8851 X 10-5 


1.1299 X 10-^ 
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1.443 X 10-5 


1.8105 X 10-5 


9.7890 X 10-5 


1.2671 X 10-^ 
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5.2846 X 10^ 


4.5133 X 10-6 


2.6891 X 10 ^^ 


3.1518 X 10-5 


0.05 


2.8857 X 10-6 


1.1439 X 10-6 


2.1081 X 10-5 


8.0360 X 10-6 


0.025 


3.7396 X 10-6 


4.3677 X 10-^ 


2.7026 X 10-5 


3.2910 X 10-6 


0.0125 


3.9559 X 10-6 


9.0239 X 10-^ 


2.8532 X 10-5 


1.1751 X 10-5 


0.00625 


4.0198 X 10-6 


3.6362 X 10-6 


2.8988 X 10-5 


4.6968 X 10-5 



Table 1: Maximum absolute error of the numerical solution compared to the exact solution (jsl) for various time 
and space discretisations h and k, with e = 0.1. 

Figure [I] shows the evolution of the numerical solution of ( 11 ) (fnum) compared with the particular 
analytical solution (JSJ) (fsoi), along with the respective error plots at each time. It's clear that for this 
range of time, the numerical solution is in excellent agreement with the analytical solution. 
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Figure 1: Evolution of numerical solution and exact solution ([3j) at (a) t = 32 {^ to) &: (b) t = 300 (ss tb), and 
the absolute errors (c) & (d) at the respective times. Parameter values are h = k = 0.1 and e = 0.1. 

Finally Table Q depicts the maximum absolute error over x for varying e at the corresponding 
t ^ ta- This measure of accuracy for the scheme is utilised when analysing the error of the weakly 
nonlinear solution for particular examples in Sections [5] and |6J 



4 Exactly solvable initial conditions 

Due to the construction of our weakly nonlinear solution it is favourable to split the initial conditions 
of the IVP Q, if possible, into an 0(1) 'exactly solvable part' (from the viewpoint of the leading order 
KdV equations) and an 0(e) perturbation. In this section we consider the case when the perturbations 
are zero. Firstly we transform the KdV equations ([7]) into the standard form 

, /^+-6/+/„^ + £^,, = 0, where (13) 
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e = /3e, V = /3r?, 
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and /3 is a free parameter. We consider the class of N-soliton solutions of equations (13) (see [16 
I19j). which yields, for our original functions, 

.^log[detM^(e,-r)], fUv.T) = 12|^log[detM^(ry,r)], 



/^(C,r) = 12^1og[detM^(e,-r)], 
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3.4930 X 10"^ 


0.0125 
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0.00625 


2024 
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Table 2: Maximum absolute error of the numerical solution compared to the exact solution ([3|) at t w t^ for 
h = K = 0.1 and various e. 

for the N X N matrix M]\}(x,T) = {mij{x,T)) with elements 

^i ~r rCj 

where ai and ki are arbitrary parameters for i,j = 1...N and 6ij is the Kronecker delta function. 

4.1 Right-propagating initial conditions 

Let us firstly consider the case of initial conditions of the IVP Q in the form of right-propagating 
N-soliton solutions of the KdV equation, i.e. 

Q2 Q3 

f\t=o = 12^1og[detMiv(a:, 0)], ft\t=o = -12^1og[detM;v(x, 0)]. (14) 

We then have the leading order solutions in the form 

/+ = 0, f- = f-{^,T). 

The higher order terms ([8]) in the weakly nonlinear solution for this case reduce to 

/^ = 0(e,r) + v(r?,r), 

where 

4>{i,T) = \l R2{x,T)dx, i:{ri^T) = -\l R2{x,T)dx, 

and R2{x,T) = {f^)\t=o, since all other terms vanish. Therefore we have 
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2 dT ./^ 



f = -oTyr / fNi^^T)dx = -6 
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-^log[detM,v(a;,-T)] 
oxol 



x=rj 



x=i 



and consequently the weakly nonlinear solution of the IVP Q for the KdV N-soliton initial conditions 
propagating to the right, is as follows 






f = 12^^1og[detMjv(e, -T)] - 6e 



02 

log[detM7v(x, -T)] 



dxdT 



\o{e'). (15) 

x=^ 



4.2 Right- and left-propagating initial conditions 

For the case of right- and left-propagating KdV N-soliton initial conditions of the IVP Q , the initial 
conditions are given by 

/|i=o = 24^1og[detM^(x,0)], /i|i=o = 0. (16) 

The leading order solutions take the following form 

/+ = /+(r?,r), /- = /-(e,r). 

Unlike the previous case the higher order terms are more difficult to determine since /"*" 7^ 0. It is 
convenient to introduce the notation 

r = 12%(e,r), f+ = 12V^r,{v,T), where 

U{x, T) = log[detMiv(a;, -T)] and V{x, T) = log[detMAr(x, T)], 
and thus the higher order terms of the weakly nonlinear solution can be written in the form 

/I = -36 [2%(C, r)K,^(77, T) + U^^^{i, T)Vr,{v, T) + K,,,(r/, T)U^{i, T)] + <A(C, T) + V^(ry, T). 

The functions (f) and ip ([9]) are constructed from the functions (10), which for this case are in the form 



Ri{x,T) = 36-^[U,,{x,T)V,{x,T) + U,{x,T)V,,{x,T)], 

R2{x,T) = l2[U^Mx,T) + V^Mx,T)] + 36^ [U^,,,,{x,T)V^{x,T) - Uc,{x,T)V^^^{x,T)] 

ox 

It therefore follows that 

R2ix,T)dx = 12[UMx,T) + VMx,T)] + 36[U^Ux,T)V^ix,T) - U^ix,T)V^^^{x,T)], 



and since we have 
1 



'P + i^ 2 



Ri{^,T) + Ri{7],T)- f R2{x,T)dx 



= 36{%(C, r)%(e, T) + U^^ir], T)V^^{ij, T) + ^^(e, T)Vi^{i, T) + V^^^{ii, r)[/,(7?, T) 
+ \ [U^tH.T) + %(e,T) - Ur,T{ri,T) - Vr,T{r,,T)]}, 

the weakly nonlinear solution of the IVP Q for N-soliton KdV initial conditions propagating both 
left and right, can be explicitly expressed in the form 

+ \ [UMx,T) + y.T(x,r)]^I«} + 0{e'). (17) 



5 Examples: weakly nonlinear solution and direct nu- 
merical simulations 

5.1 Right-propagating 1-soliton initial conditions 

We now consider particular examples of the two general classes formulated in the previous section 
and compare them to numerical simulations. Firstly for the case of initial conditions of the IVP Q 
in the form of a right-propagating single KdV soliton solution, the leading order terms are given by 



f 



0, 



r = /r(c,r) 



12^1og(l + e 



e(5-T) 



), 



k^r 



where 6{x, T) = kx + -^T + a. From (14) the initial conditions of the IVP M^ are 



/I 



,.o = 12^1og(l + e^(-°): 



Jt\t=o 



-12^1og(l + ^^^^'°^^ 



and from (|15|) the weakly nonlinear solution is 



/ = 12 



^2 

5^1og(l + e 



e(?-T')^ 



6e 



92 



dxdT 



log(i + e^(-'-^)^ 



-1 x=ri 



+ 0{^). 



(18) 



(19) 



x=5 



Explicitly evaluating the derivatives in ( 19 ) yields 

A;2„\ a 



/ 



+ 



3/c sech 
3A;4e 



k 
2 

sech^ 



e 



T + 
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3k e sech 



k^e 



sech 



k^e 



x+ 1 



1 + 



kh 



e 



k\ 



T] + 



a 



+ sech 



V 



T] + 



a 



+ Oie' 



(20) 



To consider the error of the weakly nonlinear solution for this example, we transform (20) into the 
same form used in the numerics 

3k^e^ 

+ — -- 



t + 



a 



a 



t + 



+ 0{€ 



-sech 



kV^ 



1 + 



k^e 



t + 



a 



2 

(21) 



From (18) we choose the following initial conditions for numerical simulations 
3kh sech^ f ^^! + " 1 , /, 



/, 



i,0 



i,l 



3kh sech^ r^^ 



k) + a 



where we choose some k such that the weakly nonlinear solution is applicable. 

Figure [2] depicts the time evolution of the weakly nonlinear solution (21). Figure [3^ and |3J3 
highlight the difference in the numerical solution compared with the weakly nonlinear solution, taken 
up to leading and second order, for a particular time and choice of e. The difference between solution 
(21 ) and the numerical solution, in Fig. ^, is almost indistinguishable, whilst a considerable difference 
can be observed for the leading order solution. 

Figure [3p displays the maximum absolute error for varying e at the corresponding time t ~ ta . 
As expected the errors reduce as e — )• 0, and for the weakly nonlinear solution (21), the maximum of 
the errors remains within its expected accuracy of O(e^). Each of the errors plotted in Fig. ^ are 
greater in magnitude than the corresponding errors of the numerical scheme for the same parameters 
(but for different initial conditions), given in table (pi) of Sectional Therefore we believe these plots 
provide a true measure of the accuracy of the weakly nonlinear solution and are not accountable for 
some numerical artefact of the finite difference scheme. 
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Figure 2: Evolution of the weakly nonlinear solution for right-propagating 1-soliton initial conditions, with 
e = 0.1, k = -y= and a = 0. 
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Figure 3: Right-propagating 1-soliton initial conditions, (a) Numerical solution compared with the weakly 
nonlinear solution /2 and /i, and (b) the respective absolute errors, both for e = 0.1 and t = 32 (~ ta)- (c) 
Maximum absolute errors e^^ and ej^ versus e at t ~ ta- All other parameters are k = —7= and a = 0. 
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5.2 Right- and left-propagating 1-soliton initial conditions 

The case of initial conditions of the IVP Q in the form of a right- and left-propagating single KdV 
soliton solution, yields the leading order solutions 



52 



r = /i+(r/,r) = 12^1og(l + e^('''^)), /- = /-(e,r) = 12^1og(l + e^«'-^)). 

From (16) this corresponds to the initial conditions of the IVP Q in the form 

f\t=o = 24|^log(l + e^(-'O))^ ^^1^^^ ^ o_ (22) 



The weakly nonlinear solution is therefore in the form (17) where in this particular case U{x,T) = 
log(l -|- e^(^'~'^') and V{x, T) = log(l -|- e^'^''^'). Explicitly evaluating each of the terms in the weakly 
nonlinear solution yields the solution in the following form 



/ 



3A: 



sech26'' 
sech^ei^- 



+ sedi^e''+ 



+ 



9k^e 



1 



2ax~\ 



[sech^e^ 



sech 



2ax-l^=n 



sech'^9'^+ + sech^e^^ 



- ( sech^e^- + sech^e''- 



g L ix=i 

- tanh e^- f tanh 9"^+ + tanh 0^^ 



tanh 6'^+ (tanh 9^- + tanh 9"^- 



+ 0(e'),(23) 



where we use the notation 9 = 2^{x, itT). 

We next investigate the error of the weakly nonlinear solution by firstly transforming the variables 



in (23) to the same form used in the numerics, and for numerical simulations we use the following 



initial conditions to coincide with (22) 
' y/ekx + a 



lift = 6A:^e sech^ 



fi,i = 3k'e 



sech^ (^^ 



k) + a 



+ sech2(^^ 



X + k) + a 



We again choose some suitable k to ensure the validity of the weakly nonlinear solution. 

Figure [4] depicts the evolution over time of the weakly nonlinear solution in the same form used in 
the numerics. Figure [5^ and (5)3 display the behaviour of the numerical solution compared with the 
weakly nonlinear solution, up to leading and second order, for fixed i and e. Figure [Sp displays the 
maximum of the absolute errors for the weakly nonlinear solution taken up to each order, for various 
e and the corresponding t ^ ta- Similar to the previous example the weakly nonlinear solution derived 
for this example is as accurate as expected (0(e^)), and the accuracy improves as e — )• 0. Again it's 
clear that there is a significant reduction in the maximum absolute error when the second order terms 
are included within the weakly nonlinear solution. 



11 



/ O.ID 




Figure 4: Evolution of the wealcly nonlinear solution for right- and left-propagating 1-soliton initial conditions, 
with e = 0.1, k = -y= and a = 0. 
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Figure 5: Right- and left-propagating 1-soliton initial conditions, (a) Numerical solution compared with the 
weakly nonlinear solution /2 and /i, and (b) the respective absolute errors, both for e = 0.1 and t = 32 (~ ta)- 
(c) Maximum absolute errors e^^ and e^^ versus e at i ^ ta- All other parameters are k = y= and a = 0. 
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5.3 Right-propagating 2-soliton initial conditions 



We next examine another particular example of the weakly nonlinear solution derived in Section 4.1 



namely initial conditions of the IVP Q in the form of a right-propagating 2-soliton solution of the 
KdV equation. Therefore we have the following leading order solutions 



f 



52 



'de 



r = f^{^,T) = 123^1og 1 + e^i(«'-^) + e^2(«'-^) + Ce^i(«.-^)+e2(«.-^) , 



fc?, 



where C = [{ki — k2)/{ki + k2)] and 6i{x,T) = kiX+^T+ai for i = 1,2. From (14) the corresponding 
initial conditions of the IVP Q are in the form 

/|,=o = 12^1og (1 + e^^(^'O) + e^^(^'°) + Ce^i(--0)+^^(--0)^ ^ ^^|^^^ ^ -^/l*=0' (24) 



and from (15) the weakly nonlinear solution for this example can be explicitly expressed as 



{ki - k2f + \fC (kl cosh 4" + kl cosh e\ 




D + Vc(kf cosh 0^" + kl cosh ^f 



cosh 



+ VC cosh ( ^ 



>e3 



where D = {ki — k2)'^{kf + /c|) and we use the notation 6^ 



OLi 



(where we shift 
cxi — InVC) for i = 1,2. We next examine the error of (25) by writing it in the form used for 



fc3 

kiX ± -^T + cti 



the numerics and impose the following initial conditions for numerical simulations 
{kl - /C2)2 + ^/C (kl cosh ef + k^ cosh Of^ 



fifl = 6e- 



cosh 



2 



+ vC cosh 



2 



/j,l — li,0\x=x- 



where 9f^ = kiXy/e+ai for i = 1,2, and we choose some ki and ^2 appropriately, to ensure applicability 
of the weakly nonlinear solution. 

In Fig. [6] we consider the case where the amplitudes of the leading order parts of the solution 
are close (fci = 0.61, /c2 = 0.56) and the phase shifts di = 02 = 0, such that their relative positions 
are initially close together. The initial conditions of the IVP (Wl) for this example are shown in Fig. 
[6^ and [6]i, and the evolution of the weakly nonlinear solution, up to first and second order, and the 
numerical solution, are displayed in Fig. |6J3 and [6}:, all for times within the region of validity of the 
weakly nonlinear solution. The corresponding absolute errors up to each order are displayed in Fig. |6^ 
and[6]F at each of the respective times. It's clear that the accuracy of the weakly nonlinear solution is 
within the expected accuracy of O(e^) throughout the time interval considered. There is a noticeable 
improvement in the coincidence of the numerical solution with the weakly nonlinear solution taken 
up to second order compared to the leading order solution, best emulated at earlier time (i.e. in Fig. 
[6f)and|6};). 

We also examined each combination for when the two leading order components of the weakly 
nonlinear solution are initially close or separate and when the amplitudes are similar or essentially 
different, controlled by the dj's and the kiS respectively. There were slight alterations to these results 
from the example presented in Fig. p^ however the main features remained, that being ef remained 
within O(e^) throughout the same time interval and e^ — t- e\ as time increased beyond the region of 
the validity of the weakly nonlinear solution. 



(25) 
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Figure 6: Right-propagating 2-soliton initial conditions (a) &: (d). Evolution of the numerical solution compared 
with the weakly nonlinear solution /2 and /i at (b) t = 32 (ss ta), (c) i = 300 (w tb), and the absolute errors at 
each of the respective times (e) & (f). All other parameters are e = 0.1, ki = 0.61, k2 = 0.56 and ai = 02 = 0. 



5.4 Right- and left-propagating 2-soliton initial conditions 

Lastly in this section we consider another particular example of exactly solvable initial conditions, 
namely initial conditions of the IVP Ml) in the form of right and left-propagating 2-soliton solutions 
of the KdV equations. Therefore the leading order solutions are given by 



r 



/2~(e,r) 



12|llog (1 + e^^(^'^) + e^^(^'^) + Ce^i('''^)+^^('''^)) , 
12^1og (1 + e^^«'-^) + e^2(«'-^) + Ce^i(€.-^)+e2(«'-^)] , 



and the initial conditions (16) for this example are in the form 



/|,=o = 24^1og (1 + e^^(^''0) + e^^(-'O) + Ce^i(-.o)+^2(x,o)', f^^^^^ 



0. 



(26) 



(27) 



The explicit form of the weakly nonlinear solution for this example is omitted here since the solution 



is rather lengthy, particularly due to the presence of third order derivatives of the log terms in ( 26 ) 



However the solution can be easily obtained from (17) using any computer algebra package (we used 



Maple 14)- Indeed, the only operations required in determining (17) are differentiations. Thus, 
once again we compare the weakly nonlinear solution with the numerical solution. Transforming the 



14 



variables into the same form used in the numerics we use the foUowing initial conditions for numerical 



simulations to coincide with (27) 



{ki - k2f + VC (kl cosh Of + kl cosh Of 
fi,o = 12e — : , ,„„ ,^„, , ,„„ .„„..? ^ and /j,i 



cosh 



?f'-6ifO 



+ V C* cosh 



}fO+ef 



(/. 



i,0\x=x- 



-K + /i,0|a;= 



where we again choose the /cj's appropriately in order to maintain the applicability of the weakly 
nonlinear solution. 

For the example presented here we choose the same en's and ki's as chosen in the previous example. 
The initial conditions of the IVP Q for this example are shown in Fig. uk and^, and the evolution 
of the weakly nonlinear solution, taken up to first and second order, along with the numerical solution, 
is shown in Fig. ^jp and [T]:. As the solution evolves one can notice qualitatively the resemblance of 
the leading order terms in the previous example propagating in both directions. The absolute errors 
of the weakly nonlinear solution up to each order are shown in Fig. [7^ and[7]f at each of the respective 
times and it's again clear that they are within the derived expected accuracy throughout the time 
interval considered. There is a significant improvement in the accuracy of /2 compared with /i for 
this example especially at earlier time, shown in Fig. ue, where e^/ej is approximately 0(e). 
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Figure 7: Right- and left-propagating 2-soliton initial conditions (a) & (d). Evolution of the numerical solution 
compared with the weakly nonlinear solution /2 and /i at (b) t = 32 (ps ta), (c) t = 300 (^ %), and the absolute 
errors at each of the respective times (e) & (f). Other parameters are e = 0.1, k\ = 0.61, ^2 = 0.56 and oti = 0.2 = 0. 
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6 Perturbations of exactly solvable initial conditions 

Finally we consider the weakly nonlinear solution for perturbations to the exactly solvable initial 
conditions of the IVP Q considered in Section Hi More specifically we examine just one particular 
case of initial conditions for a perturbation of a right propagating KdV N-soliton solution. The initial 
conditions of the IVP Q are chosen in the simplest perturbed form 



^2 



^3 



f\t=o = 12^1og[detMjv(x,0)] + eF^(x), ft\t=o = -12^1og[detM^(x,0)]. 

Therefore the leading order terms are still given in the form considered in Section [4. 1| 

/+ = o, /- = /-(e,r), 

but the higher order terms are of the form 

Ki{li,T) + Ki{ri,T) - 

Ri{x,T) = fHx), i?2(x,r) = ^/^(x,r). 

Therefore the weakly nonlinear solution for a perturbation of the KdV N-soliton initial conditions 
propagating to the right, is as follows 



f 



RiiCT) + Ri{7],T) - r R2{x,T)da 



where 



/ = 12^1og[detM^(e, -T)] + UfHO + FH^) " 12 



d^ 



dxdT 



log[detMiv(a;, -T)] 



x=ri 



x=e 



0(6^). (28) 



For example, we consider a perturbation of the single KdV soliton solution as in Section 5.1 and 
choose the perturbation to be defined by 

1 f kx + a^ 

F [X) = sech I — — 

Therefore the initial conditions of the IVP are 

/|,=o = 12|^log(l + e^(^'°)) + e sech ( '^^) , Mt=o 



-12£log(l + ^^^^'°^^ 



(29) 



and from ( 28 ) the weakly nonlinear solution is explicitly given by 

-7 1-1 !-!--<; .sech I — 



/ 



Sfe^sech^ i^[( 



-T] +- 
2/2 



el , ( kC + a 
H — i sech 



sech 



kr] + a 



+ 






-sech 



C - Y^] + f 1 + ^ech' 



k 



V 



—T] +- 
2 / 2 



+ 0(6^ 



(30) 



To analyse the error we transform the variables to those used in the numerics and for simulations we 

s in order 

fi^Q = Sk'^e sech^ ( " ^ "'^ ' " ) -|- e^sech 

fi,l 



use the following initial conditions in order to comply with (29) 
Skhsech'C^^^'' 






3kh sech^ ( ^^("- 



-!- — <; sech 



/ k^/ex + a 

V 2 

2 '' ^ k^/e(x — k) -|- a 



-|- sech 



k^t{x -)- k) -|- a 



2 / V 2 

The initial conditions of the IVP Q, for this example, are shown in Fig. Isk and [Sli, and the 
evolution of the weakly nonlinear solution, up to leading and second order, along with the numerical 
solution are shown in Fig. [8J3 andlSb. Figures ISfe andlSF depict the corresponding errors of the weakly 



nonlinear solution at each of the respective times. The weakly nonlinear solution (30) is within its 



derived accuracy throughout the time interval considered and is significantly more accurate than the 
leading order solution, most evident for early time (Fig. ^). However as time increases beyond the 



region of validity of the weakly nonlinear solution we find that Cj — t- e\ 
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Figure 8: Perturbation of right-propagating 1-soliton initial conditions (a) & (d). Evolution of the numerical 
solution compared with the weakly nonlinear solution /2 and /i at (b) t = 32 (ps t^), (c) i= 300 (w t^), and the 
absolute errors at each of the respective times (e) & (f). All other parameters are e = 0.1, k = y= and a = 0. 

7 Conclusion 

In this paper we have constructed a nonsecular weakly nonlinear solution of the IVP for the Boussinesq 
equation on the infinite line, for initial data generating sufficiently rapidly decaying right- and left- 
propagating waves. Seeking asymptotic multiple-scale expansions and implementing an averaging 
procedure with respect to fast time we derived two KdV equations describing the leading order terms 
and obtained formulae for the higher order corrections in an explicit and simple form. The initial data 
was split into 0(1) and 0(e) parts, and it was shown, in a case study, that this allows one to obtain 
explicit approximate solutions of the IVP for 'exactly solvable initial conditions' (from the viewpoint 
of the Inverse Scattering Transform for the leading order KdV equations). 

A finite difference scheme was implemented primarily to measure the accuracy of the developed 
weakly nonlinear solution. We showed through comparisons of the numerical solution with an exact 
solution of the Boussinesq equation that the scheme's accuracy was well within the required accuracy 
to measure the weakly nonlinear solution. Albeit that changing the initial conditions for simulations 
would undoubtedly alter the magnitude of the error of the numerical scheme determined in the test 
case in Section |3] it remains a good indicator of the scheme's accuracy. 

On comparison of the weakly nonlinear solution with relevant numerical simulations there arose two 
consistent features in each of the examples considered. Firstly, the weakly nonlinear solution remained 
within its required accuracy throughout the time interval of its validity. For a more accurate solution 
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and applicability for larger time one must reformulate equation (II| to consider higher order terms. 
Secondly, the maximum absolute error for our weakly nonlinear solution (e^) was significantly lower 
than for the leading order solution (ej). This was particularly evident for earlier times. As time 
extended beyond the end of the validity of the weakly nonlinear solution it can be seen that e^ and 
e] become comparable. 

To conclude, in all examples considered in this paper the constructed weakly nonlinear solution 
was in excellent agreement with the results of direct numerical simulations within the range of its 
asymptotic validity. It would be interesting to consider other classes of initial conditions for the leading 
order KdV equations, using the well-developed techniques of the Inverse Scattering Transform. It is 
also desirable to implement a similar procedure for water waves, which can be done along the same 
lines. Finally, one can explore the possibilities of renormalising the constructed solution, similar to 
|17j . which will be reported elsewhere. 
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